function [chi,F,C,R] = FA_em(X,r);
%%% function [F,C,R] = FA_em(X,r);
%%% estimate the factor model 
%%% x_t = C*F_t + e_t; R = E(e_t e_t')
%%% using the em algorithm 
%%% input
%%% x: TxN matrix of observables
%%% r = number of factors
%%% 
%%% output
%%% F: Txr matrix of estimated factors
%%% C: N,r matrix of estimated factor loadings
%%% R: estimated covariance matrix of the idiosincratic component e_t


OPTS.disp = 0;

[T,N] = size(X);

Mx = mean(X);
Wx = (std(X));

%%% center data
x = (X-kron(ones(T,1),Mx))*diag(1./Wx);



S = cov(x);

%%% initailize assuming the variance of the idio is equal to .5 for all
%%% variables
[v,d] = eigs(S,r,'lm',OPTS);
sigma = mean(diag(S-v*d*v'));
C = v*sqrt(d-sigma*eye(r));
R = ones(N,1)*sigma;
F = x*inv(diag(R)+C*C')*C;


err = 10;

j = 1;

while j<300 & err> 1e-4;
    F_prev = F;
    j = j+1;
    %%% E_step
    gamma = diag(1./R)*C*inv(eye(r)+C'*diag(1./R)*C);
    delta = inv(eye(r)+C'*diag(1./R)*C);
    %%% M_step
    C = (inv(gamma'*S*gamma+delta)*gamma'*S)';
    R = diag(S - S*gamma*inv(gamma'*S*gamma+delta)*gamma'*S);
    
    %%%% compute the expecep value of the factors given the estimated
    %%%% parameters
    F = x*inv(diag(R)+C*C')*C;
    err = subspace(F,F_prev); %% check convergence
end;


chi = F*C'*diag(Wx) + kron(ones(size(x,1),1),Mx); %% compute the common component

C = diag(Wx)*C; %% reattributre variance

%R = diag(diag(Wx)*diag(R)*diag(Wx)); %% reattribute variance
R1 = diag(diag(Wx)*diag(R)*diag(Wx)); %% reattribute variance
R = eye(size(R1,1));
for i = 1:size(R1,1);
    R(i,i) = R1(i);
end
